classdef GetRFMagnetometerSensitivity < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetRFMagnetometerSensitivity(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetRFMagnetometerSensitivity');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('drivingAmplitude', 0.1, @(x) x>0);         
            obj.paramParser.addParameter('signalAmplitude',0.01, @(x) x>0);
            obj.paramParser.addParameter('signalFrequency',20, @(x) x>0);
            obj.paramParser.addParameter('average', 100, @(x) x>0);           
            obj.paramParser.addParameter('order', 4, @(x) ismember(x, 1:8));
            obj.paramParser.addParameter('bandwidth', 400, @(x) x>0);
            obj.paramParser.addParameter('description', 'T=83.8degree.', @ischar);            
            obj.paramParser.addParameter('nCol',             1024, @(x) x>0);      



        end
        %%
        function obj = preprocess(obj)
            param = obj.dataPackage.parameters;
            nmrg = obj.boss;
            
            obj.boss.setRFMagnetometer('field',            nmrg.SystemInfo.MainField.bz0, ...
                'frequency',        nmrg.SystemInfo.MainField.main_frequency, ...
                'amplitude',        param.drivingAmplitude, ...               
                'order',            param.order, ...
                'bandwidth',        param.bandwidth);
            
            nmrg.bz0.setAC('frequency', param.signalFrequency, ...
                           'amplitude', param.signalAmplitude, ... 
                           'offset', nmrg.SystemInfo.MainField.bz0, ... 
                           'bandwidth', 1e3);
        end
        
        %%
        function obj = dataAcquisition(obj)
            param = obj.dataPackage.parameters;
            
            measurement_data.start = datetime('now');
            obj.boss.rbSignal.unsubscribe_signal();

            
            obj.boss.rbSignal.subscribe_signal(ziDemodSpectrum.Theta_fft_abs_avg);                                  
            measurement_data.spectrum = obj.boss.rbSignal.getDemodSpectrum(param.nCol, ...
                                                                    'average', param.average, ...
                                                                    'fSpan', param.bandwidth);                 %   [rad]  linear specturm
                                                                

            measurement_data.finish = datetime('now');
            
            obj.dataPackage.measurement_data = measurement_data;
            
        end

        %%
        function obj = analyze(obj)
            measuredData = obj.dataPackage.measurement_data;
            signal = measuredData.spectrum;
            nmrg = obj.boss;

            %------------------------
            V2nT = nmrg.bz0.setting.volt2curr * nmrg.bz0.setting.curr2nT;    % control voltage [V] to actual field [nT]    [nT/V]
            slope = 180/pi/nmrg.SystemInfo.MainField.width;                % control voltage [V] to signal amplitude [degree]   [degree/V]
            coeff = abs(slope/V2nT);                                        % [degree/nT] signal amplitude generated by 1nT transverse field in x direction

            enbw = signal.getENBW(nmrg.rbSignal.subscribed_signal{1,[1,2]});
            data = signal.getNodeData(1);
            psd = data.y /pi*180/ sqrt(enbw) /coeff * 1e3;                       % pT/sqrt(Hz)
            [ncount, edges] = histcounts(psd, 1:1:100000); % count sensitivity fro 10 ~ 100000 pT/sqrt(Hz)
            [~, loc, wid ] = findpeaks(ncount, 'NPeaks', 1, 'SortStr', 'descend');

            %------------------------
            analysis_result.enbw  = enbw;
            analysis_result.V2nT  = V2nT;
            analysis_result.slope = slope;
            analysis_result.coeff = coeff;
            analysis_result.meanSensitivity = mean( psd( psd<edges(ceil(loc+3*wid)) ) );
            analysis_result.fList = data.x;
            analysis_result.psd   = psd;
                        
            obj.dataPackage.analysis_result = analysis_result;
        end
        
        %%
        function obj = postprocess(obj)
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;
            nmrg = obj.boss;
            
            % Post operations
            nmrg.bz0.disable();
            nmrg.rbSignal.unsubscribe_signal();
            
            % summary
            summary.start = datestr(measuredData.start);
            summary.signalamplitude = param.signalAmplitude * analysis.V2nT;
            summary.signalfrequency = param.signalFrequency;
            summary.psdBackground = analysis.meanSensitivity;            
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            nmrg.SystemInfo.MagnetometorSensitivity = summary;
            
            % Store Data
            getSession;
            sess.addData('GetRFMagnetometerSensitivity', obj.dataPackage, param.description);
        end

        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            if nargin == 1
                f = figure('Name', sprintf('Sensitivity[%s]',datestr(measuredData.start)), 'WindowStyle', 'docked');
                ax = subplot(1, 1, 1, 'Parent', f);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);
                ax = subplot(1, 1, 1, 'Parent', viewer.gui.uiHandle.panelPlot);
            end
            
            fList = analysis.fList;
            psd   = analysis.psd;
            enbw  = analysis.enbw;
            sig_amp  = summary.signalamplitude;
            sig_freq = summary.signalfrequency;
            meanSensitivity = analysis.meanSensitivity;
            
            semilogy(ax, fList, psd, 'b.-', [fList(1), fList(end)], [meanSensitivity, meanSensitivity], 'r--');
            xlabel(ax, 'frequency (Hz)'); ylabel(ax, 'sensitivity pT/sqrt(Hz)');
            xlim(ax, minmax(fList(:)'));
             title(ax, sprintf('Signal = %3.2f nT @ %2.1f Hz, ENBW = %3.2f Hz', sig_amp, sig_freq, enbw));
            legend(ax, {'PSD', sprintf('Background = %3.2f pT/\\surd{Hz}', meanSensitivity)}, 'Location', 'best', 'Interpreter', 'tex');
            grid(ax, 'on');
        end

    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                startTime = arrayfun(@(t) t.start, res.summary, 'UniformOutput', false);
                if isempty(p.Results.viewer)
                    f = figure('Name', sprintf('RFResonance'));
                    ax = subplot(1, 1, 1, 'Parent', f);
                else
                    ax = subplot(1, 1, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                
                fList = cell2mat(arrayfun(@(d) d.fList', res.analysis, 'UniformOutput', false));
                psd = cell2mat(arrayfun(@(d) d.psd', res.analysis, 'UniformOutput', false));
                
                semilogy(ax, fList, psd, '.-');
                xlabel(ax, 'frequency (Hz)'); ylabel(ax, 'sensitivity pT/sqrt(Hz)');
                xlim(ax, minmax(fList(:)'));
                legend(ax, startTime, 'Location', 'bestoutside', 'Interpreter', 'tex');
                grid(ax, 'on');
            end
        end
    end
end